Clipping and intersection functions can be used to extract trajectory segments that are located within an area of interest polygon.
import pandas as pd
import geopandas as gpd
import movingpandas as mpd
import shapely as shp
import hvplot.pandas
from geopandas import GeoDataFrame, read_file
from shapely.geometry import Point, LineString, Polygon
from datetime import datetime, timedelta
from holoviews import opts
import warnings
warnings.filterwarnings('ignore')
opts.defaults(opts.Overlay(active_tools=['wheel_zoom'], frame_width=500, frame_height=400))
mpd.show_versions()
MovingPandas 0.16.0 SYSTEM INFO ----------- python : 3.10.8 | packaged by conda-forge | (main, Nov 22 2022, 08:16:53) [MSC v.1929 64 bit (AMD64)] executable : H:\miniconda3\envs\mpd-ex\python.exe machine : Windows-10-10.0.19045-SP0 GEOS, GDAL, PROJ INFO --------------------- GEOS : None GEOS lib : None GDAL : 3.5.0 GDAL data dir: None PROJ : 9.0.0 PROJ data dir: H:\miniconda3\pkgs\proj-9.0.0-h1cfcee9_1\Library\share\proj PYTHON DEPENDENCIES ------------------- geopandas : 0.13.0 pandas : 2.0.1 fiona : 1.8.21 numpy : 1.23.5 shapely : 1.8.2 rtree : 1.0.1 pyproj : 3.3.1 matplotlib : 3.7.1 mapclassify: 2.4.3 geopy : 2.3.0 holoviews : 1.14.9 hvplot : 0.8.3 geoviews : 1.9.6 stonesoup : 0.1b12
gdf = read_file('../data/geolife_small.gpkg')
traj_collection = mpd.TrajectoryCollection(gdf, 'trajectory_id', t='t')
xmin, xmax, ymin, ymax = 116.365035,116.3702945,39.904675,39.907728
polygon = Polygon([(xmin,ymin), (xmin,ymax), (xmax,ymax), (xmax,ymin), (xmin,ymin)])
polygon_gdf = GeoDataFrame(pd.DataFrame([{'geometry':polygon, 'id':1}]), crs=31256)
my_traj = traj_collection.trajectories[2]
intersections = my_traj.clip(polygon)
print("Found {} intersections".format(len(intersections)))
Found 1 intersections
ax = my_traj.plot()
polygon_gdf.plot(ax=ax, color='lightgray')
intersections.plot(ax=ax, color='red', linewidth=5)
<Axes: >
Alternatively, using TrajectoryCollection:
clipped = traj_collection.clip(polygon)
clipped
TrajectoryCollection with 2 trajectories
clipped.plot()
<Axes: >
polygon_feature = {
"geometry": polygon,
"properties": {'field1': 'abc'}
}
my_traj = traj_collection.trajectories[2]
intersections = my_traj.intersection(polygon_feature)
intersections
TrajectoryCollection with 1 trajectories
intersections.to_point_gdf()
| id | sequence | trajectory_id | tracker | geometry | intersecting_field1 | |
|---|---|---|---|---|---|---|
| t | ||||||
| 2009-02-04 10:43:04.222205 | 3225 | 773 | 3 | 2 | POINT (116.36799 39.90467) | abc |
| 2009-02-04 10:43:05.000000 | 3226 | 774 | 3 | 2 | POINT (116.36798 39.90471) | abc |
| 2009-02-04 10:43:07.000000 | 3227 | 775 | 3 | 2 | POINT (116.36795 39.90480) | abc |
| 2009-02-04 10:43:09.000000 | 3228 | 776 | 3 | 2 | POINT (116.36793 39.90487) | abc |
| 2009-02-04 10:43:11.000000 | 3229 | 777 | 3 | 2 | POINT (116.36794 39.90495) | abc |
| ... | ... | ... | ... | ... | ... | ... |
| 2009-02-04 10:48:14.000000 | 3390 | 938 | 3 | 2 | POINT (116.36535 39.90589) | abc |
| 2009-02-04 10:48:15.000000 | 3391 | 939 | 3 | 2 | POINT (116.36525 39.90589) | abc |
| 2009-02-04 10:48:16.000000 | 3392 | 940 | 3 | 2 | POINT (116.36516 39.90589) | abc |
| 2009-02-04 10:48:17.000000 | 3393 | 941 | 3 | 2 | POINT (116.36507 39.90589) | abc |
| 2009-02-04 10:48:17.399995 | 3393 | 941 | 3 | 2 | POINT (116.36504 39.90589) | abc |
170 rows × 6 columns